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Abstract 


Enhancements to an advanced ascent guidance 
algorithm for rocket-powered launch vehicles are 
described. A general method has been developed for 
conveniently and efficiently handling the common case 
of (asymmetric) launch vehicles with unbalanced thrust 
and aerodynamic moments. The new part of this 
development concerns the treatment of endo-atmosperic 
flight. An alternative method for handling the 
transversality conditions has been developed that 
eliminates the need for a priori elimination of the 
constant multipliers that adjoin the terminal state 
constraints to the performance index. As a result, new 
constraints can be formulated and implemented with 
relative ease. The problem of bum-coast-bum 
trajectory optimization is treated using a modified 
multiple shooting technique. 
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Symbols 

aerodynamic axial force magnitude 
coefficient of axial force 
coefficient of normal force 
gravity acceleration vector 
Hamiltonian function 
altitude 
Mach number 
vehicle mass 


^Guidance and flftfcfgagnitude 

flight Mechanics ^6^° P ressure 
qa product of dynamic pressure and 

angle of attack 


qP 

product of dynamic pressure and 
angle of sideslip 

R 

great circle range 

r 

vehicle position vector 

Ty, r z 

x, y, and z components of position 

s 

aerodynamic reference area 

T 

thrust magnitude 

T x 

thrust component along x-body axis 

t 2 

thrust component along z-body axis 

tco 

predicted cutoff time from guidance 

V 

vehicle velocity vector 

V*, V y , V 2 

x, y, and z components of velocity 

X 

column vector containing position 
and velocity vectors 

yb> ^b 

x-, y-, and z-body axes unit vectors 

a 

vehicle angle of attack 

p 

vehicle angle of sideslip 

r 

vehicle flight path angle 

8 

angle between velocity costate vector 
and x b vector 

<t> 

angle between position vector and 
velocity costate vector 

<t>0 

angle between Earth-relative velocity 
vector and velocity costate vector 

AT 

actual thrust minus vacuum thrust 

X 

column vector containing position 
and velocity costates 

K 

position costate 

Ky. Kz 

components of position vector costate 

K 

velocity costate 

^dir 

variable set to +1 or -1, specifying 
heads-up or heads -down flight, resp. 


parameter associated with invocation 
of linear gravity field assumption 

Subscripts 

I ,/ 

vector or quantity relative to the 
inertial frame 

r,r 

vector or quantity relative to the 
Earth-fixed frame 

0 

initial 
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Introduction 

Guidance algorithms for future reusable launch 
vehicles must be very adaptive and robust and they must 
fly near-optimal trajectories in order to meet 
challenging cost, reliability, and safety goals. For 
launch to a circular low-Earth orbit, the mission is 
performed using either a direct insertion (single burn) 
into orbit or by using a bum-coast-bum sequence.. 
Theoretically, executing two bums, i.e., the bum-coast- 
bum sequence, is more efficient, with respect to fuel 
usage, than is the single bum direct insertion. One way 
of executing the bum-coast-bum sequence is to a priori 
compute the fully-optimized trajectory and extract the 
optimal first-bum cutoff targets and the coast duration 
from the optimal trajectory. These can then be loaded 
onto the onboard computer so that the onboard 
guidance need be concerned with optimizing only a 
portion (one bum) of the overall trajectory at any given 
time. This works well and is optimal under the 
assumption of perfect modeling of the vehicle and 
environment. An increase in vehicle autonomy and 
performance can be attained by incorporating bum- 
coast-bum optimization capability into the onboard 
guidance. 

Brown, et. al. 1 , and Cohen and Brown 2 did some of 
the earliest research on multi-bum trajectory 

optimization applied to real-time guidance, using 
optimal control theory. Jezewski 3 extended the 
previous work using a linear gravity field approximation 
to reduce the computation time required to optimize 
multi-bum trajectories. Hardtla, 4 developed the 
(production) multi-bum optimization ‘Gamma 

guidance’ system for the 1US spacecraft. Unlike the 
algorithms described above, Gamma guidance used the 
impulsive bum assumption. Gath and Calise 5 
developed multi-bum optimization capability into a 
hybrid ascent guidance algorithm. 

In the atmospheric guidance research of Refs. 67 89 it 
was assumed in the formulations that the thrust was 
always directed along the x-body axis which is 
approximately true for symmetric vehicles like the 
Saturn V but not for asymmetric vehicles like the space 
shuttle. The more general case of time-varying thrust 
direction in the x-z body axes plane, applicable to 
asymmetric vehicles, is treated here. 

In all the guidance research using principles of 
optimal control theory l > 2 3.5,6,7,8,^ transversality 
conditions were handled by analytically eliminating the 
constant multipliers that adjoin the terminal state 
constraints to the cost. Each distinct set of terminal 
state constraints (each corresponding to a particular 
mission objective, e.g., circularization, cutoff at 
specified altitude and speed, cutoff at specified 
inclination and semi-major axis, etc.) requires 


independent analyses (in general) to eliminate the 
multipliers. A motivation for (analytically) eliminating 
the multipliers a priori is that it probably results in 
slightly less computation than numerical elimination. 
The motive for numerical elimination is that new 
terminal state constraint sets can be formulated with 
relative ease. A numerical approach is presented here 
that works well and incurs very little additional 
computational burden over the analytical method. 

The rest of this paper is organized as follows. The 
next section describes the trajectory optimization 
problem associated with ascent guidance. Next, the 
additional necessary conditions associated with the 
multi -bum optimization formulation are presented. The 
modifications for handling asymmetric vehicles are 
discussed next followed by a discussion of the 
numerical method for handling transversality 
conditions. Numerical results are then given using the 
guidance formulations. The paper is ended with 
conclusions and recommendations. 

Trajectory Optimization Problem 

The equations of motion for a thrusting rocket in 
atmospheric flight are: 

r = v 

t £| 0) 
m 

where the thrust components, Tx and 7z, and the axial 
and normal forces, A and N are given by: 

T X = T cos S p 

T 2 =Ts inS p ( 2 ) 

r = r VflC + AT(h) 

A = qSC A N = qSC N 

Some notes on guidance modeling follow. We’ve 
assumed that all the thrust is aligned along the x-body 
axis, x b . The velocity vector, v, can be taken as the 
Earth-relative velocity vector or inertial velocity vector 
depending on the context. For high-speed flight 
(typically occurring ‘outside’ the atmosphere), inertial 
velocity is used whereas for low-speed flight and sub- 
orbital missions, it is sufficient to let the velocity in 
question be Earth-relative velocity. Note that we’ve 
assumed the force along the y-body axis is small and 
hence we ignore it in the equations of motion and 
consequently in the optimization (but not in the guided 
simulations). In the proceeding developments, the 
position and velocity vectors will be expressed in the 
guidance coordinate frame illustrated in Figure 1. The 
guidance coordinate frame is an Earth-centered, right- 
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handed, inertial coordinate system with the x-axis 
aligned with the local vertical and the z-axis aligned 
along the expected downrange direction, it is re- 
defined each guidance cycle with the vehicle’s current 
latitude, longitude and an azimuth angle which 
approximates the downrange direction of travel. This is 
a convenient frame to work in because, for example, the 
initial position vector expressed in the guidance frame 
has y- and z-components equal to zero and the y 
components of position and velocity are typically near 
zero. 

The aerodynamic coefficients are modeled using 
least squares polynomial coefficients interpolated with 
cubic spline functions of Mach number, a common 
technique for reduction of aerodynamics coefficient 
data in trajectory optimization. The thrust difference 
term, AT, due to the effects of the atmosphere on thrust, 
can in general be represented by a cubic spline function. 
The density, p, is represented by a least squares curve 
fit of the standard atmosphere and matches the latter to 
within 5 percent up to 70 km. Similarly, a least squares 
3 rd -order polynomial fit is used to accurately model the 
speed of sound, a. 

At the time of engine cutoff, t COi k terminal state 
constraints, nonlinear functions of the states, are 
imposed: 

Vi = fi(f(w) v (w)) =0 - (3) 

Examples of terminal constraints include final position 
magnitude, flight path angle, semi-major axis, argument 
of perigee, inclination and longitude of the ascending 
node. 

In general, path constraints of the form: 

S(w)z o (4) 

are imposed where x is the state and u is the control. 
Examples include maximum axial and normal 
acceleration, maximum normal force, minimum throttle 
level, angle of attack, angle of sideslip, and min/max 
values of qa, qp and maximum dynamic pressure. The 
axial force and minimum throttle constraints are fairly 
straightforward to handle. Gath and Calise considered 
angle of attack and normal force constraints in Ref. 
Error! Bookmark not defined, using the mathematical 
rigor contained in Ref. 10, section 3.10. Constraints on 
qa and qfi can be handled in an analogous way because 
they too are functions of the state and control variables. 
The maximum dynamic pressure constraint, however, is 
a state variable constraint. It is well-known that state 
variable inequality constraints are difficult to treat using 
optimal control theory. One can also easily conceive of 


alternative, less rigorous methods of constraining 
dynamic pressure that may be just as effective and more 
appropriate for onboard guidance, e.g., Corvin 11 uses a 
feedback control law which provides maximum 
dynamic pressure control via throttle modulation. 

We consider now the question of where to point the 
z-body axis, z b> which is a function of the control x b but 
is not fully specified given x b . This decision is a 
function of how we want the vehicle to fly. We can 
choose to construct z b so that the vehicle flies at zero 
angle of sideslip: 

z b II x b x ( x b x v ) (5) 

or so the vehicle flies a zero degree (“heads-up) or 180 
degree (“heads -down”) bank angle trajectory: 

Z b II X b *( x b * r o) ( 6 ) 

The zero-sideslip option requires a non-zero roll 
angle and possibly excessive roll maneuvering. It is 
likely that future launch vehicles will have limited roll 
control authority. The heads -up/heads-d own option will 
inherently result in larger angles of sideslip but this can 
be attenuated fairly easily, if need be, by imposing a 
sideslip path constraint. In the proceeding, we adopt 
the heads-up/heads-down option. 

The optimization problem can be stated as follows. 
Determine the x-body axis history, x b (t), that maximizes 
the final vehicle mass (equivalent to minimizing fuel 
usage or minimizing flight time) subject to the 
equations of motion (1), the terminal constraints (3) and 
the path constraints (4) 

Costate Differential Equations 

We note that the atmospheric portion of flight 
occurs over a very small ground track, enabling the use 
of the flat-Earth approximations 12 : 

h = r x -r E 

g = [g 0 0 Of (7) 

With these, the state equations become (with explicit 
state dependencies called out): 
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( T i r x ) - 4rx’ V )K*„ - M r x . 

Vx = -go + 1 

m 

i T ( r x )~ A (r x , v)K^ -N{r x , v)z by (8) 

Vy ~~ 

m 

■ ... (^(^x)- 4^. v)Kz - M'x. v) z fa 

V Z 


The costate equations then are given by: 

i =i =0 

r z 

i ^ “4K * x b + ^/)^v ’ z b 

/v„ — — — “ ™ “ “ “ " 

dr x m 

i _ _ 1 ^wr^v • x b + N vx X y z b 

A vx T “ “ / Vx + 

ov Y m 

(9) 

• _ dH _ , AvyK • x b + • Z b 

A vy “ “T -A'ry + 

' m 

\ _ 3/f _ ^ , ^vz^v ' x b -^vz^v ' z b 

A vz -“T; ~ A rz + 

ov m 


where the subscripts h, vx, yy, and vz denote partial 
differentiation with respect to those variables. No 
known analytic solutions for the atmospheric 
state/costate system exist so we resort to a second-order 
Runge-Kutta numerical integration scheme to propagate 
the state/costate system. Ten integration steps are 
sufficient to obtain a good guidance solution. 

In the next subsection, we apply the maximum 
principle to obtain the optimal control for atmospheric 
flight. 

Optimality Condition 

Applying the maximum principle to the 
Hamiltonian results in the optimization sub-problem: 


j x f (7’-4 th -Aftb 


( 10 ) 


Note that the optimal control, x b ° (and hence, 
optimal z-body axis, z b °) lies in the plane defined by the 
(initial) position and velocity-costate vectors. Thus, 
with reference to Figure 2 (after Gath and Calise Error! 
Bookmark not defined.^ ^ optimization sub-problem (10) 

can be written simply as 

max {(f-^)cos(<j)+JVsin(j)} (11) 

S 


Note that we need to be able to evaluate a in terms 
of 8 so that we can evaluate A and N in the previous 


equation. From Figure 2, it is clear that a is some 
constant, a 0 , (that is, constant with respect to 8) minus 
& 

a ~a 0 -8 (12) 


The formula for ocq can be derived by setting 5 to 
zero and solving for a via: 


- v ■ x b = VCOS^ 0j | 

*b*(* b xr 0 ) __ cos(^)x b -? 0 

Z b ~ U dir ~TT\ “ ° dir ~T7\ 

sin(0) sinl©) 


v, =vz b =cr di] 


[vcos(^)cos(^j- vj (13) 


sinl 


a 0 = tan 


if <*dir[' 


VCOSI0 cosl 


-vj 


vsmld) Jcos 1 


In the preceding, note that <j> and <j> 0 are simple functions 
of the state and costate. 

The maximization sub-problem can be solved in 
many ways. One option is to take the derivative with 
respect to 8> set to zero and use an iterative procedure 
(e.g., Newton’s method) to get the root which 
corresponds to the optimum 8. This approach was 
found to be problematic because there are situations 
when the Hamiltonian (as a function of 8) is very flat 
and Newton’s method is very slow to converge. A more 
direct method is to do a Golden Section search 13 . The 
Golden Section algorithm is relatively inefficient but the 
function to be optimized in this case, the Hamiltonian, 
is fairly inexpensive to evaluate. 

Once the Hamiltonian is maximized, we need to 
construct the x- and z-body axes so that we can evaluate 
the state/costate differential equations. Start by 
expressing the x-body axis as a linear combination of 
costate and initial position vector 


x h = aik v +b\ r n 


(14) 


Dot the preceding with velocity costate and with 
initial position to get 

cos(£) = i y x b =a l +b i cos{<j>) (J5) 

cos(ftf + S) = r 0 • x b = a x cos(<f) + b x 

Solve the preceding two equations for a x and b\ in 
terms of the angles <j> and 8: 
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(16) 


( 21 ) 


a \ = 
W = 


cos(^)-cos(^)cos(^ + j) 
sin 2 (^) 

cos(^ + $)- cos(^)cos(j ) 
sin 2 (^) 


The z-body axis is constructed via; 


^l«=|M->o 

m m 

Hence, H(tj) > 0 and the end-point condition on the 
Hamiltonian can generally be replaced by the simpler 
(non-constraining) condition that the final costate 
magnitude be unity. 


z b =* b xy b 


_ °dir 

sin(^ + £) 


_ ®dU 

sinl 


in {</> + $) b 
(^cos(^)+ l\f -1 
a \K y ( a i cos(^) + ^ ) 

a \K t { a \ cos W + ) 


x h x(x b xr 0 ) 


(17) 


Transversalitv Conditions 


The transversality conditions consist of conditions 
on final costate, referred to in the sequel as ‘costate 
conditions’, and on the final Hamiltonian: 


j4/)=vV 

4/)=i 


(18) 


Typically, for a given mission, a good initial guess 
for the mission-elapsed time at which dynamic pressure 
becomes small is available. Up until that ‘simulated’ 
time, the atmospheric equations developed earlier are 
used in the (numerical) propagation of state/costate. At 
that point, the current propagated state/costate is handed 
off to the ‘vacuum propagator’ (after Ref Error! 
Bookmark not defined.). It is advantageous to do this, 
that is, only use the atmospheric equations when 
absolutely necessary, because the atmospheric equations 
are much more complicated than the analogous vacuum 
equations. Note that it is not critical that the predicted 
time of ‘atmospheric exit’, be especially accurate. 
If need be, the predicted dynamic pressure at t Exo can be 
monitored and increased if necessary. 

Numerical Solution Method 


where v is a column-vector of constant Lagrange 
multipliers, v|/ is the column-vector of terminal state 
constraints and is the ‘constraint gradients matrix’, 
that is, the k by 6 matrix whose k rows are the gradients 
(wrt state, x) of the terminal state constraints. The 
costate conditions are equivalent to requiring that the 
final costate vector be orthogonal to the space spanned 
by all admissible final state variations Sx. For a given 
set of k terminal state constraints, expressions for 6 -k 
admissible, or ‘transversality vectors’, are obtained a 
priori and the inner products of the final co state with 
these vectors are iteratively driven to zero (using a 
modified Newton’s method) simultaneously with the 
terminal state constraints. 


The multiple shooting method 14 is used to solve the 
resulting two-point boundary-value problem. The 
values of the six initial costates, and the engine cutoff 
time are the free variables that must be iterated upon to 
null the k terminal state constraints, the 6-k costate 
conditions and the constraint on final costate magnitude. 
Multiple shooting allows the user to guess state/costate 
values at more than just the initial point, significantly 
reducing the well-known sensitivity to the initial 
costates. In fact, it has been found that insertion of just 
one additional shooting point (besides the initial 
shooting point), placed just after peak dynamic 
pressure, dramatically reduces the sensitivity compared 
to single shooting. 


Note that, because the costates can be arbitrarily 
scaled by a positive factor, the end-point condition on 
the Hamiltonian is equivalent to requiring that the 
Hamiltonian be positive. 

hence, [ v T , g(r) T ] T is a transversality vector whenever 
all the terminal constraints are Keplerian constants 
(also, whenever true anomaly is free) and thus 
satisfaction of the costate conditions imply that the 
‘Keplerian’ part of the Hamiltonian (at the final time) is 
zero. Now, note that the sign of the non-Keplerian part 
of the Hamiltonian is of definite sign: 


Abort Guidance Formulation 


The abort guidance formulation is very similar to 
the nominal guidance formulation. In the case of abort 
to downrange or abort to launch site, the target values 
are obtained from an entry profile 15 , i.e., a table of 
reference altitude, speed, and flight path versus range- 
to~HAC ‘break points’, available from onboard data 
used (or generated) by the entry guidance function. An 
additional constraint for aborts to a landing site is that 
fuel be depleted during the ascent bum so that, e.g., 
landing gear loads are not exceeded. 
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A set of terminal constraints for a downrange abort 
case is given by: 
r-r d =0 

\-r-r d v d sin(^)=0 
v ■ ((r x I hac )x r)- v d cos(y d )r 2 sin(<9) = 0 
v-(rxr HAC ) = 0 

where the HAC radius vector, i* H ac> is a unit vector 
directed from the center of the Earth to the heading 
alignment cone, the desired values (‘d’ subscripts) are 
obtained from the re-entry profile, and the range angle, 
6 \ is the angle subtended by the vehicle position vector 
and the HAC radius vector. The first constraint fixes 
the altitude, the second and third fix the vertical and 
horizontal speeds, resp., and the fourth nulls the vehicle 
heading error with respect to the HAC. Note that the 
third ensures that the vehicle is actually headed toward 
the HAC and not away from it. 

For the case of downrange aborts, we want to bum 
all the propellant so, instead of optimizing fuel usage, 
we choose to maximize final speed v. (28) 
and 

To sum up the constraints, we have the four final 
state constraints (26), the transversality conditions (33), 
and the (non-constraining) constraint |X V ((/)| = 1. The 
parameters to be iteratively determined are the six 
initial co states and the time of flight. 


Numerical Results 


Will show plots of gimbal, attitude angles and 
multi -bum trajectories, altitudes, compare performance 
with single-bum trajectories showing benefits of multi- 
bum optimization. Will compare with POST generated 
open-loop trajectories. 


Summary and Conclusions 


This paper describes an ascent guidance algorithm 
that re-optimizes the entire ascent trajectory each 
guidance cycle from liftoff to main engine cutoff. High- 
fidelity guided trajectories compared with POST open- 
loop trajectories demonstrate that it provides near- 
optimal performance despite flat-Earth simplifications 
of the costate equations. The high-fidelity trajectory 
simulator with the guidance cycling at 1 Hz runs 
significantly faster than real-time, indicating the 
algorithm’s efficiency. An abort- to-do wnrange site 
formulation is given along with guided abort trajectory 
results. 

Future work will involve using day-of-launch winds 
in the guidance solution and incorporating Mach- 
scheduled angle of attack constraints. Several candidate 
RLVs are 2 -stage vehicles with gimbal ed engines, 


possibly requiring modification of the assumption that 
all thrust is directed along the x-body axis. 
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